AD-B128  261 
UNCLASSIFIED 


AN  ANALYSIS  OF  A  FINITE  ELEMENT  METHOD  FOR 
CONVECTION-DIFFUSION  PROBLEMS  <U>  MARVLAND  UNIV 
COLLEGE  PARK  LAB  FOR  NUMERICAL  RNALVSIS 
U  G  SZVHCZAK  ET  AL.  MAR  S3  BN-1002 


1/1 


NL 


F/G  12/1 


0 


INSTITUTE  FOR  '  ;  iT'f  -A  -TT:^ 
AND  TECI-INOLCXVy* 


& 


Laboratory  for  Numerical  Analysis 
Technical  Note  BN-1002 


AN  ANALYSIS  OF  A  FINITE  ELEMENT  METHOD  FOR  CONVECTION-DIFFUSION  PROBLEMS 
PART  II:  A  POSTERIORI  ERROR  ESTIMATES  AND  ADAPTIVITY 


W.  G.  Szymrzak 


rZ 

MAY  1  6  1&3 


I.  BabuSka 


«*■  a* 


fat  P*H***J 


83  05  16 


019 


March  1983 


UNIVERSITY  OF  MARYLAND 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dote  tutored) 


REPORT  DOCUMENTATION  PACE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  PEPOPT  NUMBCN 

Technical  Report  BN- 1002 

*.  RECIPIENT'S  CATALOO  NUMBER 

4.  TITLE  (end  Subtitle) 

AN  ANALYSIS  OF  A  FINITE  ELEMENT  METHOD  FOR 
CONVECTION-DIFFUSION  PROBLEMS.  Part  II:  A 
POSTERIORI  ERROR  ESTIMATES  AND  ADAPTIVITY 

S.  TYPE  OF  REPORT  4  PERIOD  COVERED 

final  life  of  the  contract 

T.  AUTMOPf»> 

W.  G.  Szymczak  and  I.  BabuSka 

*.  CONTRACT  OR  GRANT  HUMBERT*; 

ONR  N00014- 7 7-0623 

s.  pcnfopming  organization  name  ano  aoopcss 

Institute  for  Physical  Science  &  Technology 
University  of  Maryland 

College  Park,  MD  20742 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  4  WORK  UNIT  NUMBERS 

*1.  CONTPOLLING  OFFICE  NAME  ANO  aoopcss 

Department  of  the  Navy 

Office  of  Naval  Research 

Arlington,  VA  22217 

IS.  REPORT  OATS 

March  1983 

IS.  NUMBER  OF  PAGES 

34 

14.  MONITOPINC  AGENCY  NAME  4  AOORCSSfM  dlllerent  Iron  Controlling  Olll ct) 

IS.  SECURITY  CLASS,  (ol  tbit  report) 

ISa.  OECL  ASSl  FI  CATION/  DOWN  GRADING 
SCHEDULE 

14.  DISTRIBUTION  STATEMENT  (ol  i bit  Report) 

% 

Approved  for  public  release:  distribution  unlimited 

■  7.  DISTRIBUTION  STATEMENT  (ol  the  ebetrect  tntotod  In  Block  30,  II  dlllereni  In*  Rtptrt) 


it.  supplement  any  notes 


IS  KEY  WORDS  (Continue  on  retette  •Ido  II  ntc«l«r  ond  Identity  by  block  nimber) 


20.  ABST PACT  (Continue  an  retette  tide  II  necettery  end  Identify  by  blotk  number) 

A  posteriori  error  estimates  are  derived  for  the  finite  element  method  presented 

in  Part  I.  These  estimates  are  proven  to  have  the  property  that  the  effectivity 

index  0  *  (error  estimate/true  error)  converges  to  one  as  the  maximum  mesh  size 

goes  to  zero.  An  adaptive  mesh  refinement  strategy  is  based  on  equi libriating 

local  error  indicators  whose  sum  comprises  the  global  error  estimate.  Numerical 

results  show  that  0  is  nearly  one  even  on  coarse  meshes,  and  that  optimal  meshes 

are  created  by  the  adaptive  procedure.  The  successful  solution  of  a  non  linear 

problem-modeling  flow  through  an  expanding  duct,  makes  evident  the  robustness  of 
VKo  .  . .  . .  -  ...  ■  - 


do  1473  COITION  OF  I  NOV  SS  IS  OBSOLETE 

,  S.  N  0107- LF- 014- 6601 


f CCUPITV  CLASSIFICATION  OF  THIS  PAGE  fBkw*  Dele  tntered) 


■  ■>  * 


An  Analysis  of  a  Finite  Element  Method  for  Convection-Diffusion  Problems. 

* 

Part  II:  A  Posteriori  Error  Estimates  and  Adaptivity 
W.  G.  Szymczak+  and  I.  Babuska^ 


r: 


h: 

I 

|?\ 

I 


E. 


Laboratory  for  Numerical  Analysis 
Technical  Note  BN-1002 


^Applied  Mathematics  Branch,  Code  R44,  Naval  Surface  Weapons  Center,  White  Oak, 
Silver  Spring,  Maryland  20910. 

^Institute  for  Physical  Science  and  Technology,  University  of  Maryland,  College 
Park,  Maryland  20742. 


This  work  was  supported  in  part  by  the  Office  of  Naval  Research  under  contracts 
N00014-77-C0623  and  N00014-83-WR30147 . 


Abstract 


A  posteriori  error  estimates  are  derived  for  the  finite  element  method 
presented  in  Part  I.  These  estimates  are  proven  to  have  the  property  that 
the  effectivity  index  0  =  (error  estimate/true  error)  converges  to  one  as  the 
maximum  mesh  size  goes  to  zero.  An  adaptive  mesh  refinement  strategy  is 
based  on  equilibriating  local  error  indicators  whose  sum  comprises  the 
global  error  estimate.  Numerical  results  show  that  0  is  nearly  one  even  on 
coarse  meshes,  and  that  optimal  meshes  are  created  by  the  adaptive  procedure. 
The  successful  solution  of  a  non  linear  problem-modelling  flow  through  an 
expanding  duct,  makes  evident  the  robustness  of  the  method. 


CHAPTER  1 


INTRODUCTION 

This  paper  is  the  second  part  of  a  two  part  series  in  which  an  adaptive 
finite  element  method  for  convection  diffusion  problems  is  derived  and  analyzed. 
In  the  first  paper  an  upwinded  finite  element  method  was  described  which  was 
shown  to  yield  a  quasi-optimal  approximation  to  the  exact  solution  of  the 
model  problem 

-eu"  +  a(x)u'  +  b(x)u  =  f  in  (0,1) 

(1.1)  u(0)  =  a 

Bju’d)  +  B2u(l)  =  B. 

Although  quasi-optimality  is  mathematically  important,  the  goal  of  a 
numerical  computation  is  to  provide  an  accurate  approximation  to  the  exact 
solution  of  a  mathematical  model  describing  some  physical  phenomena. 

Therefore,  a  numerical  computation  should  not  only  produce  an  approximate  ‘ 
solution,  but  also  an  estimation  of  its  accuracy. 

The  effectiveness  of  such  an  error  estimate  can  be  measured  by  the  , 

effectivity  index 

0  =  (error  estimate/true  error). 

Using  the  upwinded  finite  element  method  described  in  Part  I  to  solve  (1.1), 
an  a  posteriori  estimator  is  found  such  that  0  -*■  1  as  the  maximum  mesh  spacing 
h  =  h(A)  -*•  0.  An  error  estimator  with  this  property  is  called  asymptotically 
exact.  Computationally,  the  value  of  | 1  —  © )  is  shown  to  be  small,  often  less 
than  .2,  even  on  coarse  meshes. 


These  estimates  are  derived  by  first  considering  a  projection  of  the 
true  error  onto  a  space  of  functions  which  is  zero  at  the  nodes.  This  pro¬ 
jected  error  is  shown  to  be  close  to  the  true  error  in  norm.  Since  the 
projected  error  is  zero  at  the  nodes  it  can  be  approximated  locally.  There¬ 
fore,  each  interval  has  a  local  error  indicator  associated  with  it,  and  the 
global  error  estimate  is  simply  a  sum  of  these  indicators. 

The  error  indicators  lead  to  a  procedure  for  adaptive  mesh  refinement. 

It  has  been  shown  by  Babulka  and  Rheinboldt  [1] ,  [3]  that  the  equilibriation 
of  the  error  indicators  leads  to  an  optimal  mesh.  A  strategy  based  on  this 
face  has  been  developed  and  described  in  [2] ,  [4] ,  [5] .  That  strategy  is 
used  in  this  paper  as  well. 

This  refinement  strategy  is  completely  automated  by  the  computer.  First, 
an  initial  solution  is  obtained  on  an  unrefined  mesh.  Error  indicators  are 
calculated  for  each  element  and  an  error  estimate  is  computed  by  summing  the 
indicators.  If  the  estimate  is  below  some  prescribed  tolerance,  the  algorithm 
stops.  Otherwise,  a  threshold  value  is  computed  and  all  elements  having 
indicators  above  this  threshold  are  subdivided.  The  algorithm  continues  in 
this  way  until  either  the  specified  tolerance  is  attained  or  computer  resources 
are  used  up.  Numerical  results  show  this  algorithm  to  be  very  effective  in 
resolving  boundary  layers  or  other  singularities  in  the  solution. 

An  adaptive  procedure  based  on  a  posteriori  error  estimates  has  been 
developed  by  Reinhardt  ( [13] ^  [14])  for  a  norm  which  arises  from  symmetrizing 
the  bilinear  form.  The  method  of  symmetrization  was  first  studied  by  Barrett 
and  Morton  in  [  6]  and  [7].  Unfortunately,  for  problems  of  the  form  (1.1) 
with  b(x)  4  0,  the  method  described  by  Reinhardt  requires  the  solution  of  a 
full  matrix  equation,  instead  of  the  usual  band  matrix  for  the  conventional 


or  upwinded  methods. 
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In  Chapter  2  of  this  paper  we  summarize  the  results  from  Pa'rt  I.  This 
provides  the  mathematical  framework  for  the  problem.  The  a  posteriori 
estimates  are  derived  in  Chapter  3  and  are  proven  to  be  asymptotically  exact. 
Finally,  some  numerical  results  are  presented  in  Chapter  A.  The  examples 
were  selected  in  order  to  show  the  optimality  of  the  method,  the  reliability 
of  the  error  estimates,  and  the  robustness  of  the  algorithm.  The  robustness 
is  displayed  through  the  successful  solution  of  a  non-linear  problem  with  a 
turning  point.  This  problem  arises  from  a  model  of  flow  through  an  expanding 
duct . 

The  asymptotic  theory  presented  in  this  paper  is  based  on  the  assumption 
that  the  maximum  mesh  size  h  is  small.  In  the  adaptive  mode,  this  is  not 
quaranteed.  Ratios  of  the  maximum  step  size  to  the  minimum  step  size  often 
exceeded  250  in  the  numerical  experiments  of  Chapter  A.  This  suggests  that 
this  requirement  is  not  necessary  in  practice. 

Numerical  results  for  a  simple  2-D  problem  in  which  the  flow  is  in  the 
direction  of  the  x-axis  have  been  performed  (with  a  posteriori  estimates  and 
adaptivity)  and  appear  very  promising.  These  results  as  well  as  further  1-D 
examples  are  in  preparation. 


CHAPTER  2 


PRELIMINARY  RESULTS 

In  this  chapter  we  summarize  the  results  of  ^9]  which  are  needed  in 
this  paper.  The  following  notations  are  used: 

H^(I),  k  =  0,1,...,  1  P  <  00  is  the  usual  Sobolev  space  on  the 

interval  I  =  [0,1]  consisting  of  function; 
with  k  derivatives  in  Lp(I). 

A  =  {0  =  x  <  x,  ,  <  .  .  .  <  x,  =  1} ,  where  N  =  N(A)  is  an  arbitrary 

01  N 

mesh  on  I. 


h.  =  x.  -  x.  ,  and  I, 
J  J  j-1  j 


xj)  for  J 


!)*••)  N • 


(tu  +  h^+1)/2  for 


j  =  1,...,  N  -  1, 

A  variational 


c„  =  h„  and  h  =  max  h. 

N  N  j  j 

setting  will  be  presented  for  the  problem 


Lu  =  -  cu"  +  a(x)u'  f  b(x)u  =  f  in  I, 

(2.1)  u (0)  =  •„ 

i'(u)  =  B^u'Cl)  +  B2u(1)  «  B. 

For  the  operator  L  we  assume 

Al:  a(x)  C  [0,1] ,  a(x)  >_  a  >  0, 

b(x)  C  C° [0,1] ,  b(x)  >  b  , 

2 

and  b  is  such  that  a  +  4eb  =  y  >  0. 

For  the  boundary  operator  F  we  assume 

A2:  B1,S2  _>  0,  8r  +  62  >  0. 

The  assumptions  Al  and  A2  are  sufficient  to  ensure  that  a  maximum  principle 
for  (2.1)  holds,  and  the  Green’s  function  for  L  is  bounded  uniformly  bv  a  con¬ 


stant  independently  of  e. 


We  also  have  additional  assumptions  on  the  input  data,  namely. 


A3  i) 


ii) 


The  source  term  f  is  of  the  form  f  =  f  +  f, ,  where  f  £  L, (I )  and 

o  1  o  1 

N-l 

£  =  I  C.5(x  -  x.),  where  6(x  -  x. )  is  the  Dirac  delta  function 
1  i=l  1 


N-l 

at  the  meshpoint  x  =  x  .  Furthermore,  Z  [c. | 

i=l 

of  N  and  f  is  independent  of  c. 
a  is  bounded  independently  of  e.  Also,  if  8^  ^ 


=  K  is  independent 


0  then  4^  is 

Pi 


61 

bounded  independently  of  e  and  if  8  =  0  then  —  is  bounded 

1  65 


independently  of  e. 

A4  a(x)/  £  Cr+2(I.),  b(x)/:  £  Cr+1(I  ),  fQ(x)/  £  Cr+1(I  )  (r  will  be 

j  J  i  j 

specified  later),  and  a(x)  and  b(x)  are  independent  of  e. 

In  [19]  Chapter  2  it  was  shown  that  assumptions  A1-A3  are  sufficient  to  ensure 
the  existence  of  a  unique  solution  to  (2.1)  which  is  bounded  independently  of 
e.  Assumptions  A1-A4  are  assumed  to  hold  throughout  this  paper.  Furthermore, 
without  loss  of  generality,  we  restrict  ourselves  to  the  case  of  homogeneous 
essential  boundary  conditions,  i.e.  a  =  0,  and  if  8^  =  0  then  8=0.  This 
restriction  is  made  for  the  theory  only,  and  not  for  the  numerical  examples 
presented.  Because  the  Green's  function  is  bounded  it  follows  that  if  A3  holds 
then  the  solution  to  (2.1)  is  bounded  independently  of  e.  The  assumptions  A1-A4 
are  assumed  to  hold  throughout  this  paper. 

Let  L*  denote  the  formal  adjoint  operator  to  L,  i.e., 

L*  =  —  e  ~  a(x)  +  (b-a ' ) (x) . 

dx 

The  boundary  operator  adjoint  to  f  is  F*,  where  for  u  sufficiently  smooth 


\  e  +  a (1)J  u(l)  +  Eu'd), 


if  r  0, 


u(l). 


if  =  0. 


We  now  define  the  spaces  and  bilinear  form  used  to  pose  (2.1)  variationallv 
The  space  ^  1  p  °=  is  defined  as  the  completion  of 

H1  ={uChJ(I):  u(0)  =  0,  u(l)  =  0  if  62  -  0  }, 

with  respect  to  the  norm 


(2.2)  Mu  Lo 


-  ± 

| u  [  Pdx  +  £  p  |u(x, )  | 


3  3 


,  1  <  n  <  00  , 


|u| |L  (I)  , 


where  =  N  -  1  if  ^  =  0  and  Nj  =  N  if  g1  +  0. 

The  space  H°^  can  be  easily  identified  with  Lp  ft  RN1,  that  is, 

u  ’  Hp.i  ‘  LP  *  rNi’  and 


1  n  X/P 

(I)  +  *  C.lldjl  1  ’  1  £  P  <  - 


(2.3)  j  | u I | Ho 


max [ |  | u | | L  (I  ^  |d.  |]  , 


V  oil 

In  consistencv  with  our  definition,  we  sav  u  C  H  H  H  (I)  if  u  G  H  (I)  and 

P » A  p  P 

d^  =  u(xj)  f°r  J  =  1»-  •  •»  . 

The  space  ,  is  defined  bv  .  =  (vGH^(I):  v(0)  =  0,  v(l)  =  0  if  B. 

q .  £ ,  A  *  q » c  •  ( <  q  i 


and  vL  CH‘(I.),  3  =  1,...,  Nl  for  1  <  q 


For  anv  v  C  H' 

q 


define 


bv 


( 


q .  t , 


(2.4) 


|v||h2 

q,c,A 


=  < 


N  /•  N-l  , 

Z  [L*v|qdx  +  Z  tq|J(v'(x 

j  =  l  1  j=l  1  ’ 

+  hjj-q|r*(v)  |q 


l/q 


-l 


max  [  max  ! ! l*v 1 1  ,  max  c|j(  K4))!f. 

l<j<N-l  1 


1  <  j  <  N 


V 


r*(v)|h“1] .  q  = 


On  H°  .  x  , ,  where  —  +  — 
p,A  q  1 1.:  p  q 


=  1,  1  p  <  »,  we  define  a  bilinear  form 


( • , • )  by 


(2.5)  B  (u 


N  r 

•v)  *  ,E,  Ji 


N-l 


uL*vdx  -  Z  ed.J(v'(x.))  +  d  T*(v), 

•  i  >/  x  .  ,  .  1  1  N 

1=1  i  j=i  J 


where  J(v’(x.))  =  v'(x.+0)  -  v' (x.-O)  for  1  <  j  <  N  -  1,  and  v' (x.  +  0)  =  lim  v' 
J  J  J  -  -  J  “  x-x± 

j 


The  following  results  were  proven  in  fl9] 


Lemma  2.1.  (See  1 1 9 J  Lemma  2.6).  Let  v  £  YT  then 

-  q,e,A, 


sup 


M  Ih2 

q’r’A  UCH°  A 
p,A 


I  Ba(u,v) | 


|u| lH° 

p,  A 


Theorem  2.2.  (See  fl9l  Theorem  2.8).  I  f  v  £  ,,  then  vC  L  (I)  fl  H^I)  with 

—  n  C  oo  n 


(x) . 


) 


(:./) 


IV'IIL  (I)  1  V  1/q~l  |  |v|  Ih2 

q  q  ,  F  ,  A 


i  1  q 


where  and  are  independent  of  v,  q,  c  and  A. 

For  the  finite  dimensional  trial  space,  from  which  the  approximation  is 

taken,  we  use  S  =  (u  C  C°  fl  H°  :u|T  is  a  polynomial  of  degree  <  r}.  For 
r  p ,  A  I .  — 

3 

the  test  space,  two  possibilities  are  considered.  First,  consider 


SL(n)  =  Span 


(n) 


3  !»•••*  ^*2*  ’  ~  2 


where 


(2.8a)  L* 


(n) 

-1,3 


(ri) 

'  'i.r 


n  .  . (x) ,  on  I .  and  I . , , 

-1,1  3  3+1 


0, 


elsewhere , 

for  i  ,.i  =  1,  . . .  ,N. 


(n) 


( 


(2.8b)  L*v „  .  =  \ 

X  *  J 


(n) 


x-x .  , 
J-1- 
h. 

J 


0, 


+  n  (x) ,  on  I 
-  >  J  J 


elsewhere , 


,(X.)  =  0,  for  !  =  0, . .  . ,  r-2 ,  i,j  =  1,...,N. 

<■  1 3  1 


Let  n . 


max 


£=-l, . . . ,r-I 


l|ni,.i<X)MW«p, 


and  n  =  max  n. 

3  J 


The  following  result  guarentees  that  the  finite  element  approximation 

u  C  S  obtained  when  S  ^  ^  is  the  test  space  is  quasi-opt imal  in  H°  .  if  n 
r  l  p  *  d 


is  sufficiently  small. 
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Theorem  2.3  Let  n  be  sufficiently  small  and  —  +  —  =  1.  Then  there  exists  a 
-  p  q 

Dq  >  0  bounded  away  from  zero  independently  of  e,  such  that 

inf  sup  |BA(u,v)  |  _>  Dq  ,  for  1  P  £  00  • 

v  C  S  (n)  u  C  S 
L  r 


I  lv!  iH2  =  1  IHIho  =  i 

q .  t ,  a  p  ,  a 

It  is  shown  in  [19]  Chapter  4  that  the  functions  |  can  be  found 

explicitly  such  that  n.  <  Ch . .  Because  these  functions  are  of  the  form 

j  -  ; 

/  -Ay  a0h 

P^(y)  +  P2(y)e  J  with  A  =  — —  ,  when  A  is  large  special  quadrature  rules  are 
needed  to  control  the  quadrature  errors.  In  order  to  avoid  this,  the  space 
SL(n)  is  projected  onto 


(2.9)  S  v  ;  =  Span  )  y 


J  — 1 ,  ...,  ^1*  — 1 ,  ...,  r — 2 


where 


v'  =  P  ((  (n)n  -  V  t 

xm  pk((^.j  }  }  ~  :n  ai  i  • 


with  denoting  the  projection  operator  onto  the  first  k  +  1  Legendre 

polynomials  <t>  ,...,  on  each  interval  I  .  The  following  results  holds 
OK  j 

(k) 

when  is  the  test  space. 

fkl  7 

Theorem  2.4.  Let  Sv  be  defined  in  (2.9)  with  k  =  2r  +  1.  For  v  C 
-  a  q,e  ,A 

2 

and  f  G  (H  )’  satisfying  A3,  define 
q »  e  ,  a 


F  (v)  5  <f,V>  - 


0  if  =  0, 


j  ^v(l)  if  B,  *  0, 

Then  the  following  hold:  (See  [19]  Chapter  2  for  a,  and  Theorems  3.3  and  5.3 
for  b,  c,  and  d). 


a)  There  is  a  unique  solution  u  G  H°  .  to  B  (u,v)  =F  (v)  for  each 

p,  A  A 

v  C  H2  -  l<q  <»,-+-=!,  and  u  G  H°  A. 
q,e,A  pq  ,  A 


9 


b)  There  exists  an  hc,  independent  of  c  such  that  for  all  h 


there  exists  a  unique  u  C  S  satisfying 

a  r 


B.(u,v)-F(v)  for  each  v  C  S 
A  o  Cl  t  a  a 


and  a  unique  u^C  satisfying 


B  (u  ,  v,  )  =  F  ( v.  )  for  each  v.  C  S 
ALL  L  L  L 

c)  I  I  u  —  u  I  I  o  ±  C  inf  ||u-w||Ho  ,  l£p<_°°. 

p,A  w  C  Sr  p,A 

d)  -  uJ!Ho  ,  C  max  h  r+2  {||f(r+2)||L  ( 

P.A  J  *  J 


+  a 


(r+2) 


(V 


+  b 


(r+1) 


ll»  <V 


for  1  £  P  1.  °°»  and  r  1. 

(3) 

For  the  computations  performed  the  spaces  and  were  used  as 

the  trial  and  test  spaces.  With  these  spaces.  Theorem  2.4  yields 


(2.10)  |  | u  -  uj  |Ho  <_  C  inf  |  |u  -  w|  |Ho  +  max  C  h 
P.A  wCSj^  P.A  j  J  J 

where  u  is  the  piecewise  linear  finite  element  solution, 
a 

Numerical  quadrature  was  performed  by  taking  the  piecewise  cubic 

interpolant  of  a(x)  and  piecewise  quadratic  interpolants  of  b(x)  and  f(x), 

and  then  integrating  exactly.  Using  this  quadrature  rule  the  finite  element 

solution  ii^C  exists  fi  r  h  sufficiently  small,  and  the  estimate  (2.10) 

holds  with  u  replaced  by  u,  .  The  details  of  this  result  arc  omitted  here  but 
a  h 

are  proven  in  Hgj  for  the  general  case  using  as  the  trial  space  r  ^  1. 
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CHAPTER  3 


ASYMPTOTICALLY  EXACT  ERROR  ESTIMATORS 


(3,1)  S.p.A  =  {u  "  (u»d1....»dN)e  H°>a  :  di  =  0,  i  - 


(3.2)  =  {v  CH2  :  v(x  )  =  0,  j  = 

with  1  <  p,q  <  ®,  and  —  +  —  =  1.  These  spaces  are  denoted  by  the  subscript 

—  -  p  q 

B  to  represent  the  "bubble"-like  characteristic  of  functions  which  vanish  at 
each  nodal  point. 

First,  we  consider  an  approximation  to  e  =  u  -  u^,  where  u^C  is 
the  finite  element  solution  obtained  using  the  test  space  S^n^ ,  with  n 
sufficiently  small.  Denote  by  Pe  the  solution  of 


(3.3)  L(Pe)  =  L(e)  =  f  -  L(uL)  in  I  , 


Pe(xj_1)  =  Pe(Xj)  =0,  j  =  1, . . . ,N. 

Then  Pe  also  solves  the  variational  problem:  find  Pe  £  ,  such  that 

B,p,A 


(3.4)  B  (Pe,v)  =  B  (e,v)  for  each  v  £  1C  ,  1  <  q  <  °°. 

A  A  o ,  CJ ,  A 

Because  of  assumption  A4,  Pe  C  H°  ^  and  is  bounded  independently  of  e.  The 
following  result  shows  that  u^  +  Pe  is  a  superconvergent  improvement  of  u^ 

Lemma  3 . 1  Let  e  =  u  -  u^  and  Pe  be  defined  by  (3.3)  or  (3.4).  Assume 
n is  sufficiently  small  (independently  of  e).  Then 

!  !e  ‘  T-«\  L°  1  c  n|  | e  |  Lo  ,  1  <  p  <  <*>, 
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where  C  is  independent  of  e.  A,  and  n. 

Proof:  By  (2.3),  (2.4)  and  Lemma  2.1  it  follows  that 


|u||H°  =  sup 

P,A  vCH2 

q,e,A 


!BA(u,v)| 

q»c,A 


l  <  P  < 


Therefore, 


(3.5)  I |e  -  Pe| jo  =  sup 

P,A  vCH2 

q,c,A 


BA(e-Pe,v) | 

I  M  Ih2 

q,e,A 


2  /  \ 

For  a  given  v  C  H  .,  let  w  £  S,  be  such  that 
q, e , A  v  L 


N 


(3.6) 


(n) 


w  (x)  =  E  v(x.H  1  (x), 

V  j  ~l  J 


(n) 


where  .  .  is  defined  by  (2.8a).  Clearly,  w  (x.) 

V  J 


v(x^)  and  hence 


v  -  w  £K_  . .  Therefore, 

v  B,q, A 


(3.7) 

Also,  wv  £  S  ^ 

(3.8) 

Equations  (3.6) 


B.(e-Pe,v-w  )  =  0. 
A  v 

implies  that 

B  (e,w  )  =  0. 

A  v 

(3.8)  imply  that 


(3.9)  |BA(e-Pe,v)|  =  |BA(Pe,wv)| 

N  f 

-  I  E  \  (Pe)L*(w  )j  <  2n||v||  l|Pe|| 

1 j  V 


Using  (3-5)  and  (3.9)  with  Theorems  2.3  and  2.2  we  obtain 


t  I  e-Pe  |  |  Ho  <_Cn!|Pe||L  ,  1  1  P  £  “• 

p,  A  P 

The  desired  result  now  follows  provided  n  is  sufficiently  small. 

As  mentioned  earlier  n  =  0(h),  and  hence  if  E  =  J  |Pe| ,  then  E  is 

P 

an  asymptotically  exact  estimator  to  l I e I  I jj°  ^  Lemma  3.1.  However,  we  seek 

p,A 


an  estimator  for  e  =  u-u  where  u  is  the  finite  element  solution  obtained  by 

a  a  a 

(k) 

using  as  the  test  space. 

Let  Pe  denote  the  solution  to  the  problem:  find  Pe  C  k!  .  such 
a  a  B,p,  A 

that 


2 

B  (Pe  ,v)  =  B  (e  ,v) ,  for  each  v  C  IL  1  <  q  <  ®. 

A  a  A  ot  B ,  q ,  A , 

Pe^  also  solves  (3.3)  with  e  replaced  by  e^.  Hence  PeaC  H°  By  the 
triangle  inequality  we  have 


(3.10)  llea-Pea||Ho  —  ! 1 ea  "  e I t R°  +  J  |e  -  Pe ) | Ro 

P,  A  p,  A  P,  A 


w(Xj_1)  =  w(Xj )  =0,  j  =  1,...,N. 

It  follows  from  the  maximum  principle,  specifically  see  Lemma  4.1,  that 
if  hj  is  sufficiently  small,  then 
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& 


h.  2h 

(3.11)  I  lwl  IL(b(I  )  -  C  min(-^»-p)  I  |L(uL-ua)  |  lLoo(Ij)  . 

On  1 .  , 

J 

L(u--u  )  =  -e (u  -u  )"  +  a(x)  (u  -u  )'  +  b(x)(u  -u  ), 

La  La  La  La 

and  hence  by  the  local  inverse  theorem  (since  uT -u  £  S  ) ,  we  have 

L  a  r 

(3.12)  I  lL(uL~ua)  ML  (j  )  1  Clj?  +  h~  +  I  I ul~  Uot  I  I L  (I. )  * 

“  J  hj  J  00  J 

Combining  (3.11),  (3.12),  and  Theorem  2.4d  we  obtain 


(3.13) 


Pe  -  Pe  L  <  C  max  h . 
a  L  —  .  j 

“  J 


where  C  depends  on  a(x),  b(x),  and  f(x),  but  is  independent  of  e,  and  h. 

Beginning  with  (3.10),  and  applying  Theorem  2.4,  Lemma  3.1,  and 
inequality  (3.13),  we  obtain  the  following  result. 

Theorem  3.2.  Let  e  =  u  -  u  ,  where  u  is  the  finite  element  solution 

-  a  a  a 

(k) 

obtained  when  using  S  ,  with  k  =  2r  +  1,  as  the  test  space.  Let  Pe 
°  a  a 

be  the  projection  of  e  into  ,  defined  by 

J  a  B,p,£ 

Vpea,v)  =  BA(ea,v)  for  each  v  CkJ^,  1  <  q  <  -. 

Then  by  assumption  A. 3,  Pe  £  H°  ,  and 


>,A’ 


r+2 


I ea  “  ^ea  I  I h°  -C1r'lleaUH°  +  C2  max  h  ,  for  1  <  p  <  *. 

'  P,A  j 


p,  A 


If  a  lower  bound  on  the  error  of  the  type 


(3.14) 


leJIH0  21  ch 

p,A 


r+1 


14 


€, 


i 


& 


is  available,  where  C  can  be  chosen  independently  of  c  and  A,  then  it 
follows  from  Theorem  3.2  that 


1  - 


|Pe  °  |e  -Pe  | |  0 

1  a 1 'H  .  1  a  or  'H 

- iuA\  < 


leaNH°  . 

P,A 


|eaMH°  . 
p,  A 


—  <  CjTl  +€2^1. 


Since  n  £  Ch,  it  follows  that  E  =  j ) Pe  | | ^o  is  an  asymptotically  exact 

error  estimator. 

In  order  to  justify  (3.14),  it  suffices,  for  example,  to  suppose  that 

there  is  an  interval  I  C  I  for  which 

o  ~ 

uCC°(To),  but  u  C  srd0), 

and  there  exists  a  constant  C  >  0  such  that  for  all  meshes  A  in  some  class 

min{h  :  I.  C  I  }  >  Ch(A). 
j  j  -  o  - 

For  further  details  see  [12].  In  practice  these  conditions  are  not  very 
demanding,  and  no  restrictions  in  the  algorithm  need  be  imposed. 

We  must  now  determine  a  way  to  compute  | |Pe[ |  o  *  E. 

01  H  A 

p,  A 

This  value  E  need  not  be  computed  exactly  (in  general  this  is  not  possible), 
but  It  could  be  approximated  by  some  value  E  provided  that  E  is  also  an 

A  A 

asymptotically  exact  estimator. 

Let  w  *  Pe  .  Then  w  solves 
a 

-ew"  +  aw'  +  bw  =  p  on  I. 


=  w^xj)  =  °» 

where  p  ■  f  -  Lu^  is  the  residue.  If  we  rescale  this  problem  to  I  ■  [0,1], 
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using  the  notation  g(y)  =  g(x^  ^  +  h  y)  f°r  y  C  [0,1],  and  dropping  the 
index  j,  the  function  w(y)  solves 


(3.15)  L^w  =  --^-  w"  +  w'  +  bw  =  £  in  I, 
h 

w(0)  =  w(l)  =  0. 

Consider  the  case  when  r  =  1,  i.e.,  linear  trial  functions.  Let  i>A 
be  the  solution  to 


(3.16)  Lh(wA)  = 


e  ,  1 

72  WA  +TwA 
h 


=  P,  in  I, 


wA(°>  =  w^1)  =  °* 


where  =  a(l),  and  is  the  linear  interpolant  of  p.  The  function  «A(y) 


is  calculated  explicitly  in  (3.31).  Assume  that  h^^||w 


A1  'L  (I) 
P 


w 


A'  vv 


=  ij  can  also  be  calculated  exactly.  The  value  will  be  used  as  a  local  error 
indicator  and  the  estimator  EA  is  calculated  by 


E 


A 


1/P 


In  order  to  prove  that  E  is  an  asymptotically  exact  estimator  the 

A 

following  results  are  needed. 

Lemma  3 . 3  The  Green's  functions  G^(x|y)  and  G^(x|y)  for  the  operators  L^  and  L^ 
are  positive  and  satisfy  the  following  inequalities: 


(3.17)  Gh(x|y)  <  Ch(l-e-'h/E)/Y,  for  all  x,y  C  [0,1], 
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(3.18)  Gh(x|y)  <  h(l-e 


for  all  x,y  C  [0,1] , 


)/&1, 

and 

-a.h/4e 

(3.19)  G^(x|y)  >_  yh(l-e  )/4a^,  for  0  <_  y  <_  x, 

,  1  3 

and  —  <  x  <  y, 

4  —  —  4 

-2  -  1/2 

where  C  is  independent  of  e  and  h,  y  =  (a^  +  4ct))  '  ,  a^  =  a(l),  and 

t>  =  min  b(y) . 

Proof:  Inequalities  (3.17)  and  (3.18)  follow  directly  from  rescaling  the 
bound  on  the  Green  function  established  in  [19],  Pro°f  of  Theorem  2.5. 
v  in  (0,x),  G  (x|y)  satisfies 


^Ch(K|y)  5  -  -^(Gh)yy  -  ^  <ah)y  -  0  , 

with  boundary  conditions 

Gh(x|0)  =  0,  and  Gh(x|x)  =  g(x). 

where 


,  -ha. /4c 

Let  z(y)  =  [1  -  e  ].  Then 

4al 

L*z(y)  <  0 


From  the  maximum  principle  it  follows  that  Gfxly)  >  z(y)  if  7  <  x  <  7,  and 

h  —  4  —  —  4 

0  1  y  1  x. 

Theorem  3.4  Let  w  and  w^  be  as  defined  in  (3.15)  and  (3.16),  respectively. 

Assume  that  the  approximate  solution  u  is  piecewise  linear,  i.e.,  r  =  1. 

a 

Then 


w  -  wA 


II  (I)  —  ^ih ' 

P 


'Lp(l) 


C2h 


Proof:  Since  our  attention  will  always  be  focused  on  the  rescaled  functions 
on  the  interval  I  =  [0,1]  the  "tilda"  notation  is  dropped  for  this  proof.  Fi 
a  function  g(y)  on  I  we  use  the  notation  =  g(0+)  and  g^  =  g(l~),  and  gj  i< 
the  linear  interpolant  of  g,  i.e.  gj  (y)  =  g0  +  y^1“80)  •  To  avoid  confusion 
with  fQ(x)  (see  Assumption  A3)  we  write  ^(x)  =  fQ(x)  in  this  proof. 

Let  z  =  w  -  wA.  Then  z  solves 
A 


+  bz 


P 


PI  - 


(a(y)-a1) 


h -  WA  '  b(y)wA’  in  l’ 


z(0)  =  z(l)  =  0. 


Write  z  =  z.  +  z„  +  z_  where 


First,  Consider  z^(y) .  By  assumption  A4,  a(x),  b(x),  and  |((x)  are 
piecewise  smooth.  Also  u^(x)  is  piecewise  linear,  and  bounded  (Theorem  2.4). 
Therefore, 

(3.23)  (p-PjMy)  =  h2[KA(y)  +  u\cB(y)] 

where  K^(y)  and  K^(y)  are  bounded  independently  of  e  and  h.  By  (3.20)  and 
(3.23) 

(3.24)  z^(x)  =  h2  n*  (y)  +  Gh(xl  y)dy’ 

Similarly,  we  may  write 

(3.25)  rT(y)  =  pq  +  y(p1-pc)  =  pd  +  yhhA  +  u;  yb1  , 

where  y  and  y  are  independent  of  c  and  h.  From  (3.16)  and  (3.25)  we  have 

1J 


We  now  consider  two  separate  cases. 

Case  1.  1/  -  b  u  (x,  . )  I  >  4|a  u'l.  In  this  case  u'  is  bounded,  and  hence, 

o  a  j-1  1  —  2 '  o  a 1  ct 

from  (3.24)  and  Lemma  3.3, 

(3-27)  J  zA (x) |  £  Ch^  . 

Case  2.  [3  -  b  u  (x.  )j  <  -^|a  u' I  .  In  this  case 

o  oa  j-1  2  o  a 1 

P  =  1 5-  -  bu  (x.  ,)  -  a  u'l  >  -i-la  u‘  I  . 
o  1  uo  o  a'  j-1  o  a1  —2'  o  a1 

This  is  further  broken  down  into  two  subcases  where  either  | u *  I  <  1  or 

a 1  — 

|u^|  >  1.  If  |u;|  <  1  we  have  (3.27)  again.  If  |u^|  >  1,  (3.26),  (3.19), 
(3.23)  and  (3.17)  imply  that  for  h  sufficiently  small. 
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and 


| w  |!t  •  C  (l-e~alh/4C)!ir  'h 

A  L  -  ■ 

P 


(3.28) 


A'  'L 


<  C  (l-e'Yh/c)l/|h3  . 


Thus, 


'  ZA  '  I L  -Ch  IIwAHL 
P  P 


When  considering  the  term  z  ,  since  a  is  smooth, 

B 


ax  -  a(y) 


=  Ca(y)d  -  y)w; 


where  C  (y)  is  bounded  independently  of  e  and  h.  Let  t 
a 


=  1  -  y  and  define 


v(t)  by 


v(t)  =  tw^(l-t)  +  wA(l-t) 


Then  v(0)  =  0  and  v  solves 


~  v’  (t)  -  ~  v ( t )  =  q  (t)  , 
h2  h 


where  q(t)  =  -tpj(l-t) 


T  "a'1’1'  • 


Therefore, 


|v|lLi  (1-e  alh/,£)  j  j  q  |  |  and 
L1  al  L1 


(a  -a(y)  )  ,  , 

!“h -  "a^I  1  e>I|oiNLi  +  c  ll”AnLi 


From  this,  (3.21)  and  (3.18), 
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!  ZB  (x)  I  £  Ch2(l-e  aih/E)  |  If'jl  |j  +  Chl  |wA|  |L^  . 

We  must  now  show  that  either  !  I  £•  1  i  £  C  h  ^  I  |w  ]  |  /  (1-e  al  )  or  |  |  *  T  I  | .  £ 

I  Lj  A  Li  1  L\ 

This  can  be  shown  by  arguments  analogous  to  those  used  discussing  M | j  .  The 

P 

same  two  cases  are  distinguished  as  before.  In  addition  Case  1  is  split  into 

two  subcases:  Ip  I  >  2  Ip, -r  I  and  Ip  I  <  2|r,-p  |. 

1  o  —  1  1  o  1  o 1  1  1  o 

Therefore, 


(3.29)  IUBIIL  £  C1h||wA!|L^  +  C2h3  . 

From  (3.22)  and  (3.18)  we  have 

(3.30)  M2cHL  lCh||wAl|]  . 

P  1 

The  theorem  follows  from  the  fact  that  w-w  =  z.  +  zR  +  z„  and  inequalities 

A  A  B  C 

(3. 27)- (3. 30) . 

Corollary  3.5.  Suppose  that  (3.14)  holds.  Let 


where  t.  =  h^^||w  [L  . T ,  with  w  defined  in  (3.16).  Then  E.  is  an 
jjALptil  A  A 

asymptotically  exact  estimator  to  J  j e  |  |  o  for  1  £  p  £  ro,  when  r  =  1. 

ot  H  . 

P,A 


Proof :  This  result  follows  from  the  triangle  inequality  and  Theorems  3.2  and 


3.4.  Finally,  we  remark  that  Corollary  3.3  also  holds  when  the  effects  of 
quadrature  are  included.  For  details,  again  see  [18]. 
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still  remains.  From 


The  problem  of  explicitly  calculating  !lwA!lL  (j ) 


(3.16)  we  have 
(3.31) 


»A(y)  =  A 


h.a.y/e  „ 

+  Be  J  1  +  Cy2  + 


Cy  +  D  y 


where 


C  =  h. (R.  -  R  )/2 a. 
j  1  o  1 


D  =  h.R  /a,  +  (R  -  R  )/a. 
J  o  1  1  o  1 


h.a./t .  2 

A  =  -  (C  +  D) /  (l-e  J  )al 


B  =  -A 


As 


before,  a1  =  a(l)  =  a(x^),  Rq  -  R(0)  -  R(x^_j)  -  (((x..^)  “  a0uh{l.  _  b0uh^xj 


and  R^  =  R(l)  -  R(x  )  . 


When  p  =  2M,  where  M  is  a  positive  integer,  r.  =  h. | [w  | |  can  be 

3  J  A  L  (I  ) 


computed  exactly.  However,  all  of  the  numerical  results  presented  in  the  next 


chapter  are  performed  with  p  =  1.  We  do  not  know  a-priori  if  wA  will  change 


sign  or  not,  and  if  it  does,  the  zero  of  w  cannot  be  determined  explicitly. 


Therefore,  if  w  changes  sign  in  I  we  cannot  compute  | [w  | |  exactly. 
^  A 


Suppose  that  we  let 


■;  ■'.Li 


w  dv 
A  ■ 


which  can  be  easily  calculated  from  (3.31).  The  next  theorem  shows  that  with 


certain  assumptions  on  the  exact  solution  u,  if 

N 


*  * 

E  =  Z  t. 
A  .  ,  J 
3  =  1 


then  E  is  an  asymptotically  exact  estimator  to  ||eJ|  o 

A  n.  .  « 


Theorem  3.6.  Let  u  be  the  exact  solution  and  suppose  that  u"  does  not  change 


sign  on  anv  interval  1^,  j  =  Also  assume  that  u  is  sue  a  that  (3.14) 

N 

TV  II  _  I  TV 

holds  for  r  =  1.  Then,  ii 


!  . ,  i  =  1 , . . . ,N.  Alsc 

J  *  i  f1  - 

l,  if  i.  =  h .  I  w 
1  > 


*  * 
and  E,  =  E  t.,  then  E.  is  an 
A  j,!  J  A 


asymptotically  exact  estimator  to  | |e  | |„o  ,  where  e^  =  u  -  u^. 

1 ,  A 


Proof.  By  the  triangle  inequality,  and  Theorems  3.2  and  3.4  it  follows  that 


(3.32) 


e  -  w.  <  e  -  Pe  +  Pe  -wa 

1  •>  A 11  „o  —  1  1  a  i'  '  o  11  j  A 1  1  o 

hi,a  hi,a  hi,a 


DjhMeJI  +C2h^ 

"l.A 


Also,  since  w  (x.)  =0,  j  =  0,...,N 
A  J 


(3.33)  Meu-wAU 

H1  * 
L  y  L\ 


>eV-  WaHLi(I)  +  Z  lertC*j)l0j. 

j=l 


Inequality  (3.32),  and  (3.33)  imply  that 


(3.34)  E  |e„(x  )|p  <  C  h||ej| 

1=1  J  J  h: 


+  c2h- 


1,A 


Let  e*  =  e  -  e  where  e  T  is  the  linear  interpolant  of  e  .  Then, 

(  i  j  I  <1  ‘ 


(3.35) 


e  -  e* 

01  1  H° 

1,A 


n1 


12  I 

l.A  J 


ICihMeJI  „  +  C2h3. 


l.A 


This  shows  that  |  e*  |  is  an  asymptotically  exact  estimator  for  | |e  ||  o 


Since  e*  does  not  change  concavity  on  I.,  and  e*(x.  ,)  =  e*(x.)  =  0, 
a  J  a  J-i  x  J 


""vy  '  U. e* 


>1^(1)  '  ea' 


-  &  \A 


j.i  ji. "  a  ' 


.  'M 

J=1  l*/  J 


e*  -  w. 

>J>  A 


i  l|e„-wAllLi(I)il|e,  -„Al|Li+||aJILi 


i.  Cxh| |ea| |  0  +  C2h3, 

hi,a 

where  the  last  inequality  followed  from  (3.32)  and  (3.35). 


The  assumption  uM  doesn't  change  sign  locally,  if  violated,  will  in 
general  not  disrupt  the  effectivity  of  the  error  estimator  E*.  If  u"  does 
change  sign  in  some  interval  1^,  then  1^  cannot  lie  within  the  boundary  layer 
or  interior  layer  of  the  solution  u.  This  is  because  in  the  boundary  layer 
region,  |u"|  ^  Ce  1,  where  C  is  independent  of  e  (see  c.g.  [10]).  Suppose 
that  u  is  smooth  in  this  interval  and  that  u"  is  bounded  in  I  j  independently  of 
e.  Then,  since  u"  vanishes  at  some  point  in  I ^ , 


e*(x) 


=  |u(x)  -  u  (x)  I  <  Cty  |  u"  (f  )  |  <  Ch^,  for  *CIj. 

i  J  x  J  J 


is  an  asymtotically  exact  estimator  by  Theorem  3.6,  and 


The  value 


also,  by  assumption  (3.14), 


N  *  2 

£  ll«„llL  iCh2. 

3=1  1 


Since 


it  C| 

2*11^  (I.)  —  ’  t^e  error  on  this  interval  is  one  higher  order  than 


can  be  expected  by  the  best  approximation  to  u  by  a  linear  function  on  I  . 
Thus,  in  general,  the  error  in  this  interval  is  negligible  in  its  contribution 
to  the  total  error  or  an  asymptotically  exact  error  estimate.  Even  with  this 
consideration,  the  following  precautionary  measure  is  taken  just  in  case  this 
assumption  is  violated. 

Let 


Q( 


O  *  hj  j  y*  wA(t)dt  -  J  wA(t)dt 


By  (3.31)  wA  can  have  at  most  one  zero  in  the  open  interval  (0,1).  If  v.’A ( C ) 


=  0 


and  l  C  (0,1),  then  Q(0  =  hj  I  I^a^L 


=  T  . 


Also  note  that  Q ( 1 )  =  x  .  Let 


hj I  I «A I  I Q  =  max  Q(Ci),  where 

i=l * • • • *  4 


Cj  =  1/3, 


C2  =  1/2 


=  2/3,  and 


^4  =  X- 


Then  ||*||q  is  a  norm  over  the  space  of  functions  wA  defined  in  (3.31). 

Since  this  space  is  only  two  dimensional  the  norms  |  |  •  |  |  and  |  |  ■  \  |.  are 

Q  L! 

equivalent  on  this  space.  Therefore,  if  t'  =  h.  |  |w  |  |  ,  then  C  <. 

3  3  A  i 


» 

<  C„,  and  furthermore  the  constants  C,  and  C.  are  independent  of  e 
Tj  “  2  12 
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and  h.  Thus,  even  if  the  assumption  on  u"  is  violated,  we  can  still 


quarantee  that  the  estimator  E'  = 


t!  satisfies 
1 


Cl± 


e; 


1  C2 


1,A 


with  and  independent  of  c,  and  A.  For  the  computations  presented  in 

the  following  chapter,  t!  =  h.||w.||n  was  used  as  the  error  indicator  for  I.. 

J  1  A  y  J 

In  conclusion  of  this  chapter  we  state  a  superconvergence  result  for  the 

errors  at  the  nodal  points. 

Corollary  3.7:  (Nodal  superconvergence) 

The  nodal  errors  e  (x.)  satisfy 
a  j 


N 

E 

j=l 


*„<y 


l/p 


1  c1h  | 


+  c2h- 


P ,  A 


1  <  p  <  °°. 


Proof:  For  p  =  1  this  is  (3.34).  For  p  >1,  the  result  follows  by  appropriately 
modifying  (3.32)  and  (3.33). 


CHAPTER  4 
NUMERICAL  RESULTS 


As  mentioned  in  the  introduction,  the  equilibriation  of  the  error 
indicators  leads  to  an  optimal  mesh.  An  effective  algorithm  for  performing 
this  task  has  been  developed  by  Babuska  and  Rheinboldt, and  implemented  in  FEARS 
(Finite  Element  Adaptive  Research  Solver)  ({2],  [4],  [11],  [1  5] ).  This  algorithm 
makes  a  prediction  on  the  error  indicators  after  a  future  subdivision.  The 
predictions  are  calculated  by  assuming  that  the  indicator  of  an  element  will  decrease 


bv  essentially  the  same  ratio  as  the  last  time  it  was  subdivided.  The  maximum 


predicted  error  is  then  used  as  the  threshold  value  for  subdivision — all 
elements  with  indicators  above  the  threshold  are  refined.  The  algorithm 
continues  to  solve  and  then  refine  until  the  desired  accuracy  is  attained. 

Example  1  is  a  typical  linear  convection  diffusion  equation  in  which 
all  assumptions  Al-AA  hold.  A  detailed  description  of  the  effectivity  index 
0  and  the  rates  of  convergence  for  both  uniform  and  adaptive  meshes  is 
presented.  Example  2  is  a  linear  problem  with  a  turning  point  which  violates 
assumption  A1 .  Nevertheless,  this  problem  was  succesfully  solved  by  the 
algorithm.  Example  3  is  a  non-linear  turning  point  problem.  The  robustness  of 
the  method  Is  best  displayed  through  the  results  obtained  for  this  problem.  In 
all  examples  linear  elements  were  used  for  the  trial  space. 

Example  1:  Consider  the  problem: 


(A.l)  -eu"  +  u'  +  (1  +  c)u  =  -g-ea  +  (1  +  c)(a  -  g)x  in  (0,1) 


u(0)  =  u(l)  =  0 
-(1 +c)/e  -1 

where  a  =  1  +  e  and  6  =  1  +  e  .  The  exact  solution  to  this  problem 
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1 


I 

I 

s 


I 


M 


u(x) 


(l+c)(x-l)/t  -X 
e  +  e 


+  (  i  -  fcOx. 


This  problem  was  studied  by  Kellogg  and  Han  [9]  ,  with  the  use  of  enriched 
spaces.  We  remark  that  the  results  presented  for  this  problem  were  typical 
for  all  linear  problems  satisfying  A1-A4. 

The  adaptive  process  of  computing  the  error  indicators,  refining,  and 
then  resolving,  was  initiated  on  a  uniform  mesh  of  four  elements.  The  results 
of  each  solution  pass  of  the  iteration  for  e  =  .01  are  summarized  in  Table  1. 


TABLE  1 

Summary  of  results  for  the  adaptive  mesh 
refinement  procedure  used  on  (4.1)  with  t  =  .01. 


NUMBER  OF 

MAXIMUM 

RELATIVE 

EFFECTIVIT 

ELEMENTS 

NODAL  ERROR 

ERROR  IN  H?  . 

6 

4 

2.70E-4 

1.19E-1 

1.0877 

5 

2.62E-4 

5.15E-2 

1.0442 

6 

2.62E-4 

2.27E-2 

1.0204 

7 

2.62E-4 

9.53E-3 

1.0085 

8 

2.62E-4 

4.99E-3 

1.0067 

12 

2.97E-5 

1.72E-3 

1.0052 

23 

3.91E-6 

4.68E-4 

1.0034 

35 

2.05E-6 

2.11E-4 

1.0043 

55 

2.56E-7 

7.65E-5 

1.0009 

102 

1.77E-8 

2.21E-5 

1.0005 

201 

1.18E-9 

5.69E-6 

1.0003 

One  of 

the  most  important 

aspects  of  these  computations 

is  the 

ef fectivity 

index  P  of  the  error  estimate.  Recall  9  =  E/||e 

lH°  ,  where 
1,A 
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E  is  the  computed  error  estimate  and  e  =  u  -  is  the  exact  error.  The  values 
as  listed  in  Table  1  reveal  that  the  error  estimate  at  worst  is  within  8.8% 
of  the  true  error,  and  this  percentage  decreases  rapidly  as  the  mesh  is  refined. 

Graphs  of  the  effectivity  index  are  presented  in  Figures  1  and  2  as  a 
function  of  the  number  of  elements.  Figure  1  displays  the  values  of  P  using 
uniform  meshes  for  e  in  the  range  10  ^  to  1.  Again,  notice  the  improvement 

in  0  as  the  number  of  elements  increases.  The  rate  of  decrease  of  1 1—0 |  on 

-1  -4  -2  -1 

uniform  meshes  is  0(N  )  for  e  <_  10  ,  and  0(N  )  for  e  >_  10  .  Furthermore, 

-3 

the  fact  that  the  graphs  of  9  =  ©^(e.N)  are  nearly  superimposed  for  e  10 
indicates  that  0^(e,N)  converges  as  c  +  0.  The  values  0^(O+,N)  provide  an 
upper  bound  for  0A(e,N)  on  uniform  meshes  A. 

Figure  2  compares  the  ef fectivities  when  e  =  10  ^  using  both  uniform  and 
adaptively  constructed  meshes.  With  adaptive  meshes,  the  graph  of  G>A(e,N) 
is  no  longer  smooth,  but  it  still  lies  beneath  the  graph  of  0  for  uniform 
meshes.  This  behavior  was  typical  for  all  values  of  e  tested.  The  improvement 
of  the  ef fectivities  using  adaptive  meshes  was  even  more  profound  on  other 
problems  (see  Example  2  or  [IS]). 

Next,  we  examine  the  rates  of  convergence  attained  in  the  H°  norm. 

1  ,  ii 

If  N  =  N(A)  is  the  number  of  degrees  of  freedom  in  the  mesh,  we  will  assume 

that  the  error  =  | |e| |^o  has  the  form  E^  =  CN  The  rate  y  is  assumed 

1 ,  A 

to  depend  on  both  e  and  R,  i.e.,  y  =  y(e,R),  where  R  is  the  relative  error 
.  Figure  3  shows  the  graphs  of  y(e,.l),  y(e,.01)  and  y(e>.0025), 

when  uniform  meshes  are  used.  For  R  <_  .01  the  graphs  begin  with  a  value  of 
two  and  then  decrease  to  one  (or  nearly  one),  as  e  +  0.  This  transition  from 
two  to  one  is  delayed,  as  smaller  relative  errors  are  considered.  The  optimal 


-2  - 1 

rates  of  convergence  (N  if  1/N  <  c,  and  N  if  1/N  >"•  c)  for  uniform  meshes 
are  reflected  by  these  graphs. 

_2  -3  _5 

Figure  4  shows  the  graphs  of  v (c,  10  ),  y(e,  10  ),  and  y (e,  10  ),  when 

adaptively  constructed  meshes  are  employed.  Notice  that  these  graphs  always 

lie  above  or  on  the  value  Y  =  2,  which  implies  that  the  rates  of  convergence 

_2 

observed  are  always  0(N  ).  It  is  also  important  to  note  that  when  using 

adaptive  meshes,  the  rate  of  convergence  is  unharmed  as  e  -*•  0. 

The  rates  of  convergence  using  adaptive  and  uniform  meshes  for  c  =  1,  .01. 
and  .0001,  can  also  be  seen  from  the  slopes  of  the  graphs  of  the  relative  H°  ^ 
errors,  displayed  in  Figure  5.  By  extrapolating  the  graph  of  the  error  for 
e  =  .0001  using  uniform  meshes,  it  can  be  seen  that  in  order  to  obtain  the 
same  accuracy  achieved  with  an  adaptive  mesh  of  149  elements,  approximately 
17,000  uniform  elements  would  be  required.  For  smaller  values  of  e  this  effect 
is  even  more  pronounced. 

The  maximum  nodal  errors  for  this  problem  are  displayed  in  Figure  6. 

In  each  case,  using  both  uniform  and  adaptively  constructed  meshes  for  e  =  1.0, 

_3 

and  .0001,  the  superconvergent  rate  of  0(N  )  is  attained.  Notice  that  when 

e  =  .0001,  the  maximum  nodal  errors  are  smaller  for  uniform  meshes  than  for 

adaptive  meshes.  This  unexpected  result  does  not  contradict  the  theory  —  the 

adaptivity  optimizes  the  mesh  based  on  the  errors  in  without  regard  to 

1 ,  A , 

the  nodal  errors. 

To  describe  the  distribution  of  the  mesh  points,  we  use  the  mesh  grading 
function  £^(x).  The  function  £^(x)  is  the  piecewise  linear  function  on  the 
mesh  with  the  property 

=  3  —  0,  .  .  • ,N. 

For  example,  on  a  uniform  mesh,  £^(x)  =  x.  The  derivative  of  ^(x)  is  a  measure 
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of  the  density  of  the  mesh  points.  The  mesh  points  will  be  heavily  concentrated 
in  regions  where  £^(x)  increases  rapidly. 

Figure  7  shows  the  graph  of  the  mesh  grading  function  when  c  =  .01,  for 
an  adaptively  constructed  mesh  of  201  elements.  Note  the  heavy  refinement  in 
the  boundary  layer  region  (near  x  =  1).  The  remainder  of  the  nodal  points 
were  distributed  nearly  uniformly  throughout  the  rest  of  the  interval.  This 
reflects  the  fact  that  u  is  smooth  (although  not  linear)  outside  of  the 
boundary  layer. 

Example  2.  Consider  the  following  turning  point  problem: 

2 

(4.2)  -eu"  -  xu'  =  ctt  cos(ttx)  -  (irx)sin(Trx)  on  (-1,1), 

u (— 1 )  =  -2,  u(l)  =  0. 

The  solution  to  this  problem  is 

u(x)  =  cos(ttx)  +  erf  (x/*^2e)/erf  (l//2c) , 

which  has  an  interior  layer  at  x  *  0. 

Although  this  equation  violates  assumption  A1  at  x  =  0,  the  results 

presented  here  reveal  the  robustness  of • the  algorithm  and  suggest  that  the 

crucial  theoretical  results  also  hold  for  this  type  of  problem. 

Figure  8  displays  the  graphs  of  the  errors  in  the  H°  norm  for  e  =  .01 

1 ,  A 

and  .0001,  using  both  uniform  and  adaptive  meshes.  The  optimal  rate  of 
_2 

convergence  (0(N  )),  is  quickly  realized  in  all  cases.  The  errors  using 

adaptive  meshes  are  smaller  than  with  uniform  meshes,  particularly  when 
e  =  .0001. 

The  graphs  of  the  maximum  nodal  errors  for  e  =  .0001  are  presented  in 
Figure  9.  These  graphs  indicate  that  nodal  superconvergence  occurs  even  in 
the  presence  of  a  turning  point,  provided  adaptive  meshes  are  used. 

Figure  10  displays  the  graphs  of  the  effectivity  indices,  9A(e,N),  for 
e  =  .01,  and  .0001  using  both  uniform  and  adaptive  meshes.  Even  on  this 
turning  point  problem,  we  have  good  effectivity  indices— at  worst  9  =  1.72. 


Example  3:  Consider  one-dimensional  flow  in  a  duct  of  variable  cross 
sectional  area.  A  one  equation  model  for  steady  non-heat  conducting, 
viscous  flow  has  been  developed  by  Shubin  and  Stephens  in  [16]  and  [17]. 

This  model  is  governed  by  the  equation 

(4.3)  -eu  +  rx  +  G  =  0  on  (0,L) 

where 

G  =  -(y-I)C(J  -  f  )  A’/A^  + 

+  y (Du  +  E/u)A'/A^, 
r  =  J  (Du  +  E/u) , 

D  =  D  -  eA',  D  =  (y+1)C/2y,  and 
E  =  (y-1)CH/y. 

Y  =  1.4  is  the  ratio  of  specific  heats, 

C  =  .68471  and  H  =  3.5  are  constants, 

A(x)  =  1.398  +  .347  tanh  (.8  x  -  4)  is  the 
cross  sectional  area, 
e  is  the  viscosity  coefficient  assumed  to 
be  constant,  and 
u(x)  is  the  velocity. 

u 

In  the  inviscid  case  (c  =  0),  the  velocity  is  sonic  if  u  =  a  =  (E/D) 2  =  1.0801 
supersonic  if  u>a  and  subsonic  if  u<a.  When  the  viscous  problem  is  solved 
with  boundary  conditions  such  that  u(0)>a  and  u(L)<a  an  interior  layer  arises. 

In  the  limiting  case  as  e-K3,  this  layer  becomes  a  shock.  The  boundary  conditions 
used  are  u(0)  =  1.299  and  u(L)  =  .505  with  L  =  10. 

We  solve  the  nonlinear  problem  by  iterating  on  a  linearization  of  the 
equation.  This  procedure  is  described  in  [8]  and  referred  to  as  a  variant  of 
Newton-Kantorovich  method.  In  general,  consider  the  nonlinear  problem 


-eu"  +  F(x,u,u’)  =  0  in  I  =  [xo,x^] 
u(xq)  =  a,  u(xN)  =  8 

The  linearized  boundary  value  problem  for  the  Newton-Kantorovich  iteration  is 

(4.4)  -eu"  .  +  F  ,  (x,u  ,u')u'  .  +  F  (x,u  ,u')u  . 

m+ 1  u  mm  m+ 1  u  mm  m+ 1 

=  Fu(x’VUm)um+  Fu’(x’Vum)um  "  F(x’Um’Um)  in  l’ 

Um+1(xo)  =  1  * 

Vl(XN)  =  6  ’ 

The  aim  was  now  to  use  the  continuation  method  with  respect  to 

decreasing  values  of  e  together  with  intermittent  mesh  refinements  based 

on  the  linearized  equation  (4.4).  This  procedure  was  started  using  a  large 

value  of  c,  a  uniform  mesh,  and  a  linear  solution  u  which  satisfied  the 

o 

boundary  conditions.  For  sufficiently  large  z  (r  =  .1),  the  iterative  process 
converged  quickly.  Using  the  continuation  method,  the  value  of  z  would  be 
decreased  by  a  factor  of  about  two.  This  would  lead  to  convergence  provided 
the  mesh  was  sufficiently  refined  and  the  approximation  was  sufficiently 
close  to  the  exact  solution  for  the  current  c.  In  this  manner  the  adaptive 
method  not  only  found  accurate  solutions  and  error  estimates,  but  also 
significantly  increased  the  efficiency  of  the  continuations  used. 

The  computed  solutions  for  z  -  .1,  .01,  and  .001  are  graphed  in  Figure  11. 
The  exact  inviscid  solution,  which  has  a  shock  at  x  =  4.816  is  practically 
indistinguishable  from  the  computed  viscous  solution  with  e  =  .001.  The 
effectivity  indices  and  relative  errors  were  also  calculated  with  respect  to 
the  exact  viscous  solutions  and  are  shown  in  Table  2. 
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TABLE  2 


The  effective  tv  indices  and  relative  errors  for 
the  duct  flow  problem 


e 

EFFECTIVI TY  INDEX 

RELATIVE  ERROR 

.  1 

00 

. 10E-2 

.01 

.83 

. 20E-3 

.001 

.77 

. 33E-4 

The  exact  solutions  were  approximated  by  a  computed  solution  on  a  sufficiently 
fine  mesh.  The  mesh  grading  function  for  the  e  =  .001  solution  with  84  elements 
is  displayed  in  Figure  12. 

The  interior  layer  of  the  solution  contains  a  turning  point  of  the 

u 

linearized  equation  (4.4).  The  turning  point  x^  occurs  when  um(xt)  =  (E/D(xt))2 

in  which  case  F  ,(x  ,u  (x  ),u* (x.))  =  0.  If  this  point  occured  in  an  interior 
u  t  m  t  m  t 

region  of  an  element,  no  upwinding  was  performed  there.  If  F  i  was  nearly 
zero  at  one  endpoint  of  an  element  and  sufficiently  large  at  the  other  endpoint, 
then  upwinding  was  performed  based  on  the  larger  value.  This  corresponds  to 
the  "switching  schemes"  derived  in  [16]  and  [17]. 
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FIGURE  1 :  Effectivity  indices  for  problem  (4.1)  with  various  values 

of  e  using  uniform  meshes.  For  e  <_  .001,  the  values 
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FIGURE  2:  Effectivity  indices  for  problem  (4.1  )  when  e=10 

adaptive  and  uniform  meshes. 


FIGURE  3: 


Rates  of  convergence  for  problem  (4.1)  at  various  rel 
ative  errors  using  uniform  meshes. 


FIGURE  4: 


Rates  of  convergence  for  problem  (4.1)  at  various  rel' 
ative  errors  using  adaptive  meshes. 
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FIGURE  7:  Mesh  grading  function  for 


FIGURE  8: 


Errors  in  the  H?  .  norm  for  problem  (4.2) 
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FIGURE  10:  Effectivity  indices  for  problem  (4.2). 


